function [ peakAmplitudes peakTimes peakRiseTimes summationFactor response t ] = analyzeAlphaTrain( traceFile, delay, ISI, numPulses )
%ANALYZEALPHATRAIN Measures peak values of response to train of
%alpha-synapse current stimulation
%
% Measures the amplitude above baseline of each peak in the response to a
% train of alpha synapses. Assumes positive-going responses.
%
% Arguments:
%  traceFile    e.g. 'BS0090/BS0090AAAA0045.xsg'
%  delay        time to first pulse onset, in sec, e.g. 0.1
%  ISI          time between puls onsets, in sec, e.g. 0.05
%  numPulses    # of alpha-synapse input pulses
% 
% 2009-10-09 Ben Suter
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data = load(traceFile, '-mat');
sampleRate = data.header.ephys.ephys.sampleRate;
response = data.data.ephys.trace_1;

% Set up indeces
numSafetySamples = 1;
idxBaseline = 1:floor(delay*sampleRate) - numSafetySamples;
for p = 1:numPulses
    sta = ceil((delay+(p-1)*ISI)*sampleRate) + numSafetySamples;
    sto = sta + floor(ISI*sampleRate) - numSafetySamples;
    idxPulse(p, :) = sta:sto;
end

% Subtract baseline
response = response - mean(response(idxBaseline));

for p=1:numPulses
    [maxVal, maxIdx] = max(response(idxPulse(p,:)));
    peakAmplitudes(p) = maxVal;
    peakTimes(p) = idxPulse(p,1) / sampleRate;
    peakRiseTimes(p) = maxIdx / sampleRate;
end

% summactionFactor: difference in amplitude between last and first peak,
% relative to first peak. 
summationFactor = (peakAmplitudes(numPulses) - peakAmplitudes(1)) / peakAmplitudes(1);

t = 0:length(response)-1; t = t / sampleRate;
end

